%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Integrating GPE by the Method of Steepest Descent.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear;

dt = .02;
g  =100;						                         	   %Self-energy-ratio
x  =[-1:.02:1]*1.5;	   	x=[x;x];
Psi=exp(-x.^2);		    	Psi=Psi/norm(Psi(1,:));		%Trial-Wavefunction
V  =x.^2;                                             %Potential
V  =2.5*(x.^4-2*x.^2+1);

for k=1:300
  [ gradPsi,Dumm]=gradient(Psi);
  [nablaPsi,Dumm]=gradient(gradPsi);
  H=-nablaPsi + V + g*Psi.^2;
  Psi=Psi.*(1-H*dt);
  Psi=Psi/norm(Psi(1,:));

clf;			           	   figure(1);
plot(x(1,:),V(1,:),'b',x(1,:),5*Psi(1,:)+min(H(1,:)),'r',x(1,:),H(1,:),'g')
axis([min(x(1,:)) max(x(1,:)) 0 5]);
title('\fontsize{18}BEC-Wavefunction');
set(gca,'TickLength',[.03 .05],'FontSize',16,'LineWidth',1.2);
xlabel('\fontsize{18}{\itx/a_{ho}}');
ylabel('\fontsize{18}{\it\Psi/a_{ho}^{3/2}}');
pause(.01);
end;
